% by YLX on 6/10
% Version 2.3
% ~ is tilde not tilt 
%% added if structer in inner for loop 
%% use poisson_part to reduced calculate by -I*M*K
function R = findR(S, Q)
    % S is S_test
    global P J I CURRENT
    m_max = CURRENT;
    R = zeros(I+1, m_max+1); 
    % k <= m_max in the Paper. Changed to CURRENT.
    R(1, :) = 1;
    % 取 R(1,:)为初始状态 1，R(i+1,:)对应slot i，方便计算
    for i = 1:I 
        for m = 1:m_max
            for k = 1:m_max
                for k_tilde = 1:m_max
                    if k==1 % this accelerates the function
                        poisson_part = 1 - function_F(m + k_tilde);
                    else
                        poisson_part = function_Poiss(m + k_tilde - k);
                    end
                    R(i+1, k) = R(i+1, k) + poisson_part * Q(i, m) * R(i, k_tilde);
                end
            end
        end
    end
    R = R(2:end, :); % remove the first row of zeros
    return;

